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Computation  of  Hankel  Functions 
Leslie  A.    Berry- 
Methods  of  computing  solutions  to  Bessel's  equation, 

z^  — 5-  +  z  — ^  +  (z^   -  v^)  y  =  0,    for  complex  v  and  z,    and 
dz  dz 

for  small,    real  v,    are  given  in  detail.      The  methods  are 

most  suitable  for  use  on  an  electronic  digital  computer. 

With  a  computer  which  carries  eight  decimal  digits,    the 

methods  presented  for  small,    real  v  will  return  at  least 

six  correct  significant  digits  for  all  z.      When  v  and  z  are 

both  complex,    the  accuracy  depends  on    |z  |   (since  asymptotic 

forms  are  used)  and  the  ratio  v/z.      Relative  error  curves 

are  shown  as  a  function  of    |v/z  |  parametric  in    |z  |.      For 

V  sufficiently  far  from  z,    seven  correct  digits  are  obtained, 

which  for  v  —  z,    the  relative  error  is  of  the  order  10    ^  /  \z\. 


1.      Introduction 

Solutions  to  electromagnetic  or  sound-wave  problems  in  spherical 
or  cylindrical  co-ordinates  usually  involve  solutions  of  Bessel's  equa- 
tion,   z^  — X-  +  z  -r-    +  (z^   -  v^)y  =  0.      The  Hankel  functions   [Watson, 
dz  dz 

1958]  Hy     (z)  and  H^     (z)  are  one  set  of  linearly  independent  solutions. 
Other  solutions  can  be  written  in  terms  of  these  solutions   [see  (19)]. 
The  superscript  indicates  the     "kind,  "  v,    the  "order,  "  and    z,    the 
"argument"  of  the  functions,    and  either  v,    or  z,    or  both  may  be  com- 
plex.     The  purpose  of  this  paper  is  to  present  in  one  place  detailed 
instructions  for  computing  the  Hankel  functions  on  an  electronic  digital 
computer,    and  to  investigate  quantitatively  the  relative  errors  incurred. 
No  new  mathematical  development  is  involved,    but  since  considerable 
algebraic  reduction  and  systemization  of  scattered  results  were  necessary, 
publication  of  the  results  seems  desirable. 

Section  2  details  methods  of  conciputing  Hy(z)  when  v  and  z  are 
complex.     If  V  —  z,    Hi(z)  must  be  computed  accurately,    so  techniques 
are  given  in  section  3    for  computing  Hv(z)  when  v  is  real  and  small 
and  z  is  unrestricted.      The  accuracy  of  the  methods  is  discussed  in 
section  4, 


2,      Complex  Order,    v  or  z  Large 

In  this  section,    it  is  assumed  only  that  either  v  or  z  is  large. 
What  is  meant  by  "large"  will  be  discussed  quantitatively  in  section  4, 
but  the  asymptotic  expressions  used  here   [(1)  and  table  1  ]  accurately 
describe  the  qualitative  behavior  of  the  Hankel  functions  over  raost  of 
the  V  and  z  planes,    since  the  minimum  relative  error  is  about  one  per 
cent  for    |v|  —  |z|  ^l   (figure  5). 

Two  approximations  are  used.      In  the  first,   which  will  be  called 
the  Hankel  approximation  [Fok,    1934], 

H^/^{z)  s^exp  [(-1)1^-1  in/6]  Vl   -  ri  coth  T]    Hi^^(-i?),   k=l,  2,  (1) 


where 

^  =  V  (tanh  r\  -  r]),  (2) 

cosh  r\  =  v/z,    and    |arg  v|   <  tt/2   [Olver,    1954].      For    |arg  v|  >  n/2,  (3) 
(9)  is  used. 

The  principal  square  root  is  used;  i.  e.  , 

arg  vw  =  g-  arg  w, 

and  -  TT  <  arg  w  =    TT. 

for  v/z  ^1    [Fok,    1946], 

-^  =  1  0-Sr  [-10 -f^) -0-1^)1-    <^) 


The  first  term  on  the  right  has  been  used  extensively  as  an 
approximation  to  -  i§  in  radio  propagation  theory  [Fok,    1946;  Bremmer, 
1949]. 


The  second  approximation,    called  the  Debye  approximation,   will 
be  written  in  terras  of  [Watson,    1958] 


^-— ivtanhri  ^^^  ^'"       (^j  tanh  T]  J 


The  combination  of  S\     (z)  and  S^     (z)  to  be  used  depends  on  the 
location  of  v/z  in  the  complex  plane.      The  first  four  A„   are   [British 
Association  Mathematical  Tables,    1952]: 

Ao   =  1, 

15  ,2 

^^  =  8    -  24    ^°'^     '^' 

3  77  ^,2  385  ^,  4 

^^  =  128    -  576   "°*^     ^+3456   "°*     ^' 

5  1521  ,  2  17017  ,  4  17017  ,  e 

>3  -  l024    -  256^  ^°^^     ^+  I3i240    ^°*     ^  "  248832  ^°^^     ^-  ^^^ 


Forms  for  A4  --  A7  can  also  be  found  in  the  above  reference.     In 
(5),    use  the  principal  square  root.      The  series  in  (5)  can  be  written 


00 

Y    r(m+|)  A„ A^coth  11       3  Ag  coth^  ri       1  5  A3  coth^  T] 

L        r(^)      (±ivtanhTl)»  "  V  ^  v^ 


v^ 


m  =0 


(7) 

This  is  an  asymptotic  series  and  should  be  truncated  at  the  smallest 
term;  i.  e.  ,    all  terms  preceding  the  smallest  term  should  be  summed. 

The  regions  of  the  plane  to  be  considered  are   shown  in  figure  1, 
and  the  expressions  for  Hy(z)  are  given  in  table  1.      These  expressions 
are  correct  only  if  Re(z)  ^  0,    but  if  Re(z)   <  0,    the  formulas   [Watson, 
1958] 


and 


sinTTv  smTTv 


where  k  is  any  integer,    can  be  used  to  evaluate  the  desired  function. 
Although  expressions  can  be  derived  for  Hv(z)  in  regions  6b,    7b,    and 
9,    it  is  simpler  to  compute  H_y(z)   [so  that  -v/z  is  in  6a,    7a,    or  8] 
and  use   [Watson,    1958] 

h1^^(z)  =  exp  [(-1)^  vni]  h1.V(z).    k  =  1,    2.  (9) 


2.  1.     Determination  of  T] 

Directions  for  computing  T]  and  for  determining  the  correct  region 
of  the  plane,    figure  1,   will  now  be  given.     Let  T\  =  a  +i3,   where  a  and 
3  are  real.      Then  a  +  i3  =  cosh"''"  v/z.      The  infinitely  many-valued 
function,    cosh"'',    is  made  single  valued  by  requiring  that  0  ^  3  <  TT 
[Watson,    1958].      3  =  0  if  v/z  is  real  and  v/z  ^  1.     Specifically 


a=-log(|^±y^    -1    I).  (10) 


The  sign  is  chosen  so  that  a  has  the  same  sign  as  Im(v/z).     If 
Im(v/z)  =  0,    and    |v/z  j   >  1,    the  positive  sign  is  chosen;  if    |v/z  |   ^  1, 
11  =  Cos"^  v/z.      Then 


3  =  Cos"^    [Re(v/z)/cosh  a]  (U) 

These  directions  are  exact  for  use  in  (5),    and  are  valid  for  (1)  over 
most  of  the  plane.      (For  a  discussion  of  the  region  of  validity  for  the 
Hankel  approximation,    consult  Olver  [1954].)    For  (1),   when  v  =  z, 
arg  (-i§)  can  be  determined  by  examining  the  first  term  ^f  (4);  when  v 
is  not  near  z,    the  Debye  approximation  is  superior. 

The  m  in  table  1  is  the  largest  integer  less  than 
{{1   -  a  tanh  a)  tan  3  +   3}/TT. 


Now,    let 

f(a,  P)  =  1    -   e  cot  3  -  a  tanh  a  .  (12) 

The  signs  of  f(a,  0),    f(a,  0)  +  n  cot  0,    Re(v/z),    and  Im(v/z)  uniquely 
determine  in  which  of  the  regions  shown  in  figure  1  v/z  is  located. 
Coliamn  2  of  table  1   shows  how  the  region  can  be  determined. 

3.     Real  Order 

To  use  (1),    a  routine  for  computing  Hi(z)  must  be  available.     In 
this  section,   methods  are  given  for  computing  Hv(z)  for  small,    real  v. 
The  restriction  that  v  be  real  is  probably  not  necessary  if  z^  is  care- 
fully defined,    but  it  is  retained  in  this  section, 

3.  1.     Small    |z  | 
If  V  is  not  an  integer  [Watson,    1958] 

h1^^(z)  =  {-1)^  [exp  [(-1)'^  vHi]  J^(z)   -  J_^(z)3/i  sin  vTT,    k  =  1,    2, 

(13) 

where  Jv(z)  is  the  Bessel  function  of  the  first  kind,    defined  by  [Watson, 
1958] 


■v(z)  =  y 


(-l)''(z/2)^-"^'' 


L       m  !  r  (v+  m+  1) 

m  =0 

(v/2)"  r  (z/2)^  (z/2)^  {z/2)^ 

r(v+l)    V         1  !  (v+1)       2  !  (v+l)(v+2)      3  !  (v+l)(v+2)(v+3) 

(14) 

Jy(z)  is  made  single-valued  by  specifying  arg  (z/2)^  =  v  arg  (z/2)  for 
A  -  2tt    <  arg  (z/2)  ^    A.      Ordinarily,    A  =  TT,    but  may  have  different 
values  for  some  applications. 

If  V  is  an  integer,    (13)  is  indeterminate.      To  compute  Hn(z)  for 
integral  n  and  small  z,    compute  Jo(z)  and  Ji(z)  using  (14)  and  then  find 
the  Bessel  function  of  the  second  kind,    Yo(z),    from  [Watson,    1958] 


Y,(.)=|[<Y+log(f))J.(z)+f     '^(^'•^■,'^'  ].  (15) 


m=i 


where  y  is  Euler's  constant,    y  =  .  5772156649  •  •  •  .      The  Jn(z),    m  >  1 
are  found  with  the  recurrence  formula  [Watson,    1958], 


J.+  i(2)  =^  J„(z)  -  J„-i(z).  (16) 

Then  [Watson,    1958] 

y,(z)  =  {j,(z)  Yo(z)  +:^]  /Jo(z)  (17) 


and 


Finally, 


Yn,.-i(z)  =^  Y„(z)  -  Y„-,(z).  (18) 


H^^(z)  .  J„(z)  +  (-l)^-n  Yn(z).  (19) 


3.  2.      Large    |z  | 

As    |z  I  gets  large,    (14)  converges  slowly  and  loses  precision. 
More  serious  is  the  precision  lost  in  (13)  and  (19),    since  in  half  of  the 
z -plane  the  terms  on  the  right  are  nearly  negatives  of  each  other  when 
Im(z)  is  large.      For  example,    four  significant  digits  are  lost  in  (19) 
when  n  =  1,    z  =  5  exp  (i  85°),      It  is  therefore  necessary  to  use  the 
asymptotic  series   [Watson,    1958]- 

jj(l)|^j  ^    exp(iz-i(v+i)T./2)  T^,.2i^)_  ,.„<  ^^g  ,  <  2„)_         (20a) 

y  n  z/2 

jj(^3)|^j  __   exp(-iz-H(v+|)n£2)  ^^j2i.),  (-2n<argz<n),         (20b) 

^J  n  z/2 


where 


T,(z)  =  wy'^-'-^'"^-'-^''---'^-'-'^"-"''.  (21) 

^  n!  (4zr 


n  =  l 


Unless  V  =  n  +  -g,    for  integral  n,    Ty(z)  is  an  asynaptotic  series  and 
should  be  stopped  at  the  smallest  term.     If  v  =  n  +  ^,    Tv(z)  terminates, 
and  (20)  is  the  finite  series  for  Hn  +  i(z). 

Using  a  computer  which  carries  eight  decimal  digits  in  floating 
point  operations,    about  four  correct  significant  digits  in  Hy(z)  are 
obtained  over  all  of  the  z-plane,    using  the  methods  presented  above. 
This  accuracy  is  increased  to  at  least  six  digits   if  Dingle's   [1958, 
1959  ]  convergence  factor  is  used  in  (21),    as  described  below. 


3.  3.      The  Convergence  Factor 

In  a  series  of  papers,    R.    B,    Dingle  has  developed  a  general 
theory  [Dingle,    1958]  of  convergence  factors  and  has  applied  this  theory 
to  many  special  functions   [Dingle,    1958,    1959].      When  the  product  of  a 
convergence  factor,    Cb(z),    and  the  nth  term  of  an  asymptotic  series  is 
added  to  the  previous  n  -  1  terms,    the  sum  is  the  exact  sum  of  the 
entire  series.     In  practice,    of  course,    the  sum  is  not  "exact,"  but  the 
use  of  a  convergence  factor  does  significantly  extend  the  range  of 
application  of  an  asymptotic  series. 

In  the  particxolar  case  of  the  Hankel  functions,    (21)  is  written: 

N-l 

(4v^  -  l^)(4v- 

n!  (4z)" 


Y      (4v"-l^)(4v"-3")---(4v^-(2n-l)^) 
^y^^)  =  ^  +  i  n!(4zr 


nz:l 


(4vS.l3)...(4v^_(2N-l)^) 
+  N!  (4z)N ^n(v,  z),  (22) 


where  N  is  chosen  so  that  the  term  corresponding  to  N-l  is  the  smallest 
term  of  (21).      Then,    [  Dingle,     1959]  for  s  =  N  -  v  -  |, 


r  f^  .\  ~\   r(v  +  |)r(N+v-t  +  |)        t  .(t),  . 

t=o 


'     ^    '       (N+v-^)      '     ^    '       (N+v-i)(N+v-|-)      '     ^    ' 


where   [Dingle,    1958] 


(23) 


Ai°^z) 


r{s  +  l)      J      1+x/z     ' 

0 


2  A     ^z)  .  (s  +  z+l)  A^   '(z)   -  z, 


.(2) 


(i)/.x    ,    a(°), 


2zA^,   '(z)  =  (s  +  z)  A^,   '(z)  +  A^^   >{z)   -  1, 


and 


tzAl*)(z)  =  (s  +  z  +  2-t)  a[*"')(z)  +  A['~^h^),    t>2. 


(24) 


When   |s  I   ^  1,    or    |z  |   ^  1,    and    |arg  z  |   <  3TT/4  [Dingle,    1958], 


a'/N^)-     ^ 


s  +  z 


1  - 


z(s  -  2z)       z (s^  -  8sz  +  6z^) 


(s  +  z)*^  (s  +  z) 


(s  +  z)' 


•)•   (25) 


(0) 

Otherwise,    A^   '(z)  can  be   found  by  numerical  integration.      For  the 
general  routine  used  with  this  paper,    A^   '(z)  was  found  with  an  eight 
point  Gaussian  quadrature  formula: 


x^  e-'^   f(x)  dx  ^  ^    W,  f(x,), 


(26) 


^=1 


since  s  ^  5  for  use  in  (24).      The  W^   and  x^   are  taken  from  Rabinowitz 
and  Weiss   [1959].      Then  Cfg(v,  z)  is  approximated  by  the  first  four 
terms  of  (23). 
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4,      Discussion  of  Results 
Fortran  subroutines  were  written  to  compute  Hj^(z)  and  Ha(z)  and 

3  3 

Hn(z)  using  the  methods  of  section    3.      The  values  for  Hj^(z)  and  H2(z) 

3,  3 

were  checked  against  the  Tables  of  Modified  Hankel  Functions  of  Order 
One -Third  and  Their  Derivatives   [The  Staff  of  the  Computation 
Laboratory,    1945].      The  computed  values  agree  with  the  tabular  values 
to  at  least  six  significant  digits  where  the  tables  give  six  digits.      (The 
tables  give  eight  decimal  places.  )     The  values  of  Hn(z)  computed  by  the 
routine  were  checked  against  NBS  tables  for  Jo(z)  and  J3^(z),    [Mathe- 
matical Tables  Project,    1947]  and  Yo(z)  and  Y   (z)   [Computation 
Laboratory,    1950]  and  gave  comparable  results. 

The  Debye  and  Hankel  approximations  were  programmed  sepa- 
rately and  checked  for  real  order  by  comparing  the  values  they  returned 
with  the  correct  (to  at  least  six  significant  digits)  answers  given  by  the 
methods  of  section  3.     A  subroutine  valid  for  small,    complex  v  and  z 
was  obtained  from  SHARE   [Goldstein,   Kresge',    and  Chen,    I960],    and 
the  asymptotic  routines  described  here  returned  answers  correct  to 
the  expected  accuracy,    providing  an  independent  check  of  the  routines' 
logic. 

Defining 

„    ,      .  I  answer  -  approximation  i  ,_^, 

Relative  error  =       ^^ ,  (27) 

answer 

figures  2,    3,    and  4  show  relative  error  curves  for  both  approximations 
as  functions  of    Iv/z|,    parametric  in    |z  |  for    |arg  z|   =  0°,    5    ,    and  15°, 
respectively.     It  is  apparent  from  these  curves  that  the  Debye  approxi- 
mation is  the  better  except  when  v  —  z.     Even  when   |v/z  |   =  1,    if 
|arg  v/z  I  and    |z  |  are  not  small,    the  Debye  approximation  is  superior, 
as  shown  in  figure  5.      Nevertheless,    any  general  routine  for  computing 
Hv(z)  must  include  the  Hankel  approximation.      The  relative -error 
curves  all  show  that,    for  any  z,    there  is  an  €■  >  0  such  that,    if    jv-z  |   <  e, 
the  Debye  approximation  is  useless. 

If  series  (7)  is  truncated  at  the   smallest  term,    the  relative  error 
in  the  Debye  approximation  is  of  the  order  of  the  magnitude  of  that 
term.      Examination  of  figures  2  through  5  yields  an  empirical  formula. 


[relative  error  I  Sf  10    ^  I  \z\,  (28) 

for  the  Hankel  approximation  (1). 

A  general  routine  for  Hy(z)  was  written,    incorporating  both  the 
Debye  and  Hankel  approximations.      The  routine  uses  the  Debye  approxi- 
mation unless  the  smallest  term  of  series  (7)  is  larger  than  10"^/  (z  |, 
in  which  case  the  Hankel  approximation  is  used.      This  routine  and 
other  computer  routines  described  in  this  paper  are  available  from 
SHARE   [Berry,    1963], 


5.      Conclusions 

Using  an  electronic  computer,    solutions  to  Bessel's  equation  of 
real  order  and  complex  argiiment  can  be  computed  very  accurately. 
Even  when  the  order  is  complex,    computations  can  be  made  sufficiently 
accurate  for  most  applications.      The  subroutines  incorporating  the 
methods  should  be  valuable  in  radio  and  sound  wave  computations. 
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7.      List  of  Symbols 

H\     (z)  Hankel  fxinction  of  argument  z,    order  v,    and  kind  k, 

k  =  1,    2. 

Jy(z)  Be s sal  function  of  the  first  kind,    of  argixment  z  and  order 

V. 

Yy{z)  Bessel  function  of  the  second  kind,    of  argument  z  and 

order  v,    sometimes  called  Neumann's  function  and  denoted 
N,(z). 

Tv(z)  The  series  defined  in  (21). 

T]      =  a  +  ip  =  cosh"^  — 

z 

(k) 

S^     (z)  The  series  defined  in  (5). 

(k) 

Aj,  The  coefficient  of  S^     (z),    see  (6). 

f(a,  3)  =         1  -  3  cot  3  -  a  tanh  a 

m    =  Largest  interger  less  than  [(1    -  a  tanh  a)  tan  0  +   3}/Tr. 

C|s|(z)  Convergence  factor  for  asymptotic  series,    see  section  3.  3. 
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Region 

Determined 
by 

„;»,z, 

„;^',., 

1 

Ua,  p)>  0 

Ua,  P)+TT  cot  P>  0 

s;»,., 

sj^'u, 

2 

f(a.P)<0 

f(ff,  P)+TT  cot  (3>0 

Im    V  ^  <  0 

z 

s;»,.,  -  s;^',., 

s;^',., 

3 

f(«.P)<o 

f(Qr,  P)+TTCOtp>    0 

Im  .Vj  <  0 
z 

s;",., 

s;^',.,  -  s;" ,., 

4 

f(<»,  P)>  0 

f{a,  P)+Trcotp<  0 

Im  (v/z)  >  0 

S^      (z) 

„    (2)            2vTri      (1) 
S         (z)-e          S        (z) 

V                                        V 

5 

Ha.  P)>  0 

f(a,  p)+-iTcotp<  0 

Im  (v/z)<  0 

„    (1),    ,      -2vTTi„    (2),    , 
S^      (z)-e            S^       (z) 

s;^',., 

6a 

f(a,P)<  0 
f(a,P)+Tr  cotp<0 
Im  (v/z)  >  0 
Re  (v/z)  >  0 

P          (m+l)vTii     sin  mvir 
1                                    sin  VTT 

s^%)J%) 

V                  V 

S^^'(z)  +  e<"^+^>^^^      ^^"  "^^"  S<Hz) 
V                                           sin    VTT       V 

6b 

f(a,p)<0 
f(Q',  p)+Tr  cotp  <  0 
Im  (v/z)  >  0 
Re  (v/z)  <  0 

Compute  H            (z)  and 

-V 

use  (9) 

Compute  H           (z)  and 
use  (9) 

7a 

Ha,?,)  <0 
f(a,P)+Trcotp  <  0 
Im  (v/z)  <  0 
Re  (v/z)  <  0 

<i''\)    g-(m+l)vTri      sinmvTT^(2) 
V                                       sin    VTT       v 

r    ^-(m+DvTTi  sin  mvul   (2|(1) 

L                                           sin   VTT       1   V     *     '      v      *     ' 

7b 

f(a,P)  <0 

f(Q',P)+TTCOtP    <0 

Im  (v/z)  <  0 
Re  (v/z)  <  0 

Compute  H            (z)  and 
use  (9) 

Compute  H         '(z)  and 
use  (9) 

8 

Im  (v/z)    =  0 
Re  (v/z)  >  1 

1     c    (1),    >        e    (2),    V 
-  S^      (z)  -  S^       (z) 

^S<^\z).s(2\z) 
2     V                    v 

9 

Im  (v/z)  =  0 
Re  (v/z)  <  -  1 

Compute  H            (z)  and 

-V 

use  (9) 

Compute  H            (z)  and 
use  (9) 

Table  1 


Ik 


Im(v/z) 


v/z  plane 


'8— Re(V'z) 


Figure  1.     Regions  of  the  v/z  plane  (table  1). 
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Relative  errors  of  the  Debye  and  Hankel  approximations  for 
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E.elative  errors  of  the  Debye  and  Hankel  approximations  for 
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Relative  errors  of  the  Debye  and  Hankel  approximations  for 
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Figure  5.       Relative  errors  of  the  Debye  and  Hankel  approximations 
for    |v/z  1   =  1,    as  a  function  of    |arg  v/z  |,    for  various    |z 
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